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We  present  a  simulation  and  multi-objective  optimization  framework  for  the  integration  of  renewable 
generators  and  storage  devices  into  an  electrical  distribution  network.  The  framework  searches  for  the 
optimal  size  and  location  of  the  distributed  renewable  generation  units  (DG).  Uncertainties  in  renewable 
resources  availability,  components  failure  and  repair  events,  loads  and  grid  power  supply  are 
incorporated.  A  Monte  Carlo  simulation— optimal  power  flow  (MCS-OPF)  computational  model  is  used 
to  generate  scenarios  of  the  uncertain  variables  and  evaluate  the  network  electric  performance.  As  a 
response  to  the  need  of  monitoring  and  controlling  the  risk  associated  to  the  performance  of  the  optimal 
DG-integrated  network,  we  introduce  the  conditional  value-at-risk  (CVaR)  measure  into  the  framework. 
Multi-objective  optimization  (MOO)  is  done  with  respect  to  the  minimization  of  the  expectations  of  the 
global  cost  ( Cg )  and  energy  not  supplied  (ENS)  combined  with  their  respective  CVaR  values.  The  multi¬ 
objective  optimization  is  performed  by  the  fast  non-dominated  sorting  genetic  algorithm  NSGA-II.  For 
exemplification,  the  framework  is  applied  to  a  distribution  network  derived  from  the  IEEE  13  nodes 
test  feeder.  The  results  show  that  the  MOO  MCS-OPF  framework  is  effective  in  finding  an  optimal 
DG-integrated  network  considering  multiple  sources  of  uncertainties.  In  addition,  from  the  perspective 
of  decision  making,  introducing  the  CVaR  as  a  measure  of  risk  enables  the  evaluation  of  trade-offs 
between  optimal  expected  performances  and  risks. 
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1.  Introduction 

Over  the  last  decade,  the  global  energetic  situation  has  been 
receiving  a  progressively  greater  attention.  The  adverse  environ¬ 
mental  effects  of  fossil  fuels,  the  volatility  of  the  energy  market, 
the  growing  energy  demand  and  the  intensive  reliance  on  cen¬ 
tralized  bulk-power  generation  have  triggered  a  re/evolution 
towards  cleaner,  safer,  diversified  energy  sources  for  reliable  and 
sustainable  electric  power  systems  [1-6],  The  challenges  involved 
have  stimulated  both  technological  development  of  new  equip¬ 
ment  and  devices,  and  efficiency  improvements  in  design,  plan¬ 
ning,  operation  strategies  and  management  across  generation, 
transmission  and  distribution. 

In  this  paper,  we  focus  on  distribution  networks  and  the 
conceptual  and  operational  transition  they  are  facing.  Indeed,  the 
traditional  passive  operation  with  unidirectional  flow  supplied  by 
a  centralized  generation/transmission  system,  is  evolving  towards 
an  active  operational  setting  with  integration  of  distributed  gen¬ 
eration  (DG)  and  possibly  bidirectional  power  flows  [7,8], 

DG  is  defined  as  ‘an  electric  power  source  connected  directly  to  the 
distribution  network  or  on  the  customer  site  of  the  meter’  [8-10]  and 
in  principle  offers  important  technical  and  economical  benefits.  Under 
the  assumption  that  the  distribution  network  operators  have  control 
over  the  dispatching  of  the  DG  power,  improvement  of  the  reliability 
of  power  supply  and  reduction  of  the  power  losses  and  voltages  drops 
can  be  achieved.  Indeed,  DG  allocation  on  areas  close  to  the  customers 
allows  the  power  flowing  through  shorter  paths,  and  therefore, 
decreasing  the  amount  of  unsatisfied  power  demand  and  enhancing 
the  power  and  voltage  profiles.  Thus,  the  eventual  intermittence  of 
the  centralized  power  supply  can  be  smoothed  [11],  In  addition,  the 
modular  structure  of  the  DG  technologies  implies  lower  financial 
risks  [12,13]  and  thus  the  investments  on  the  power  system  can  be 
deferred  [1,3], 

Most  of  the  actual  DG  technologies  make  use  of  local  renewable 
energy  resources,  such  as  wind  power,  solar  irradiation,  hydro- 
power,  etc.,  which  makes  them  even  more  attractive  in  view  of  the 
requested  environmental  sustainability  (e.g.,  the  Kyoto  Protocol 
[7,14,15]).  Given  the  intermittent  character  of  these  energy 
sources,  their  implementation  needs  to  be  accompanied  by  effi¬ 
cient  energy  storage  technologies. 

Attentive  DG  planning  is  needed  to  seize  the  potential  advan¬ 
tages  associated  to  DG  integration,  taking  into  account  specific 
technical,  operational  and  economic  constraints,  sources  and  loads 
forecasts  and  regulations.  If  the  practice  of  selection,  sizing  and 
allocation  of  the  different  available  technologies  is  not  performed 
attentively,  the  installation  of  multiple  renewable  DG  units  could 
produce  serious  operational  complications,  in  fact,  counteracting 
the  potential  benefits.  Degradation  of  control  and  protection 
devices,  reduction  of  power  quality  and  reliability  on  the  supply, 
increment  in  the  voltage  instability  and  all  related  negative 
impacts  on  the  costs,  could  become  impediments  for  integration 
of  DG  [1-3,8,10,14,16-20], 


Viewing  DG  planning  as  a  fundamental  baseline  of  advance¬ 
ment,  many  efforts  have  been  made  to  solve  the  associated 
problem  of  DG  allocation  and  sizing.  Objective  functions  consid¬ 
ered  for  the  optimization  are  of  economic,  operational  and 
technical  type.  Among  the  first  type,  cost-based  objective  func¬ 
tions  have  been  used  considering  the  costs  of  energy  and  fuel  for 
generation,  investments,  operation  and  maintenance,  energy  pur¬ 
chase  from  the  transmission  system,  energy  losses,  emissions, 
taxes,  incentives,  incomes,  etc.  [1-3,7,8,11,13,14,16-27],  The  sec¬ 
ond  type  of  operational  objective  functions  mainly  revolves 
around  indexes  such  as  the  contingency  load  loss  index  (CLLI) 
[23],  expected  value  of  non-distributed  energy  cost  (ECOST), 
system  average  interruption  duration  index  (SAIDI),  system  aver¬ 
age  interruption  frequency  index  (SAIFI)  [7,16,28],  expected 
energy  not  supplied  (EENS)  [28,29],  among  others.  Regarding  the 
third  type  of  objective  functions,  technical  performance  indicators 
include  energy  losses  [1,30]  and  total  voltage  deviation  (TVD)  [18], 
Power  Flow  (PF)  equations  are  typically  solved  within  the 
optimization  problem  to  evaluate  the  objective  functions,  while 
respecting  constraints  and  incorporating  non-convex  and  non¬ 
linear  conditions.  Given  the  complexity  of  the  optimization  pro¬ 
blem,  heuristic  optimization  techniques  belonging  to  the  class  of 
Evolutionary  Algorithms  (EAs)  have  been  proposed  as  a  most 
effective  way  of  solution  [10],  including  particle  swarm  optimiza¬ 
tion  (PSO)  [23,24,27,31,32],  differential  evolution  (DEA)  [18]  and 
genetic  algorithms  (GA)  [3,7,11,13,14,16,26,33,34], 

An  additional  difficulty  associated  to  the  problem  is  the  proper 
modeling  of  the  uncertainties  inherent  to  the  behavior  of  primary 
renewable  energy  sources  and  the  unexpected  operating  events 
(failures  or  stoppages)  that  can  affect  the  generation  units.  These 
uncertainties  come  on  top  of  those  already  present  in  the  network, 
such  as  intermittence  and  fluctuation  in  the  main  power  supply 
due  to  unavailability  of  the  transmission  system,  overloads  and 
interruptions  of  the  power  flow  in  the  feeders,  failures  in  the 
control  and  protection  devices,  variability  in  the  power  loads  and 
energy  prices,  etc.  These  uncertainties  are  incorporated  into  the 
modeling  by  generating  a  random  set  of  scenarios  by  Monte  Carlo 
simulation  (MCS);  the  optimization  is,  then,  executed  to  obtain  the 
optimal  expected  or  cumulative  value(s)  of  the  objective  function 
(s)  under  the  set  of  scenarios  considered  [2,3,7,16,28,32,34,35]. 

In  the  search  for  the  optimal  DG-integrated  network,  the  use  of 
only  mean  or  cumulative  values  as  objective  function(s)  of  the 
optimization  hinders  the  possibility  of  controlling  the  risk  of  the 
optimal  solution(s):  the  optimal  DG-integrated  network  may  on 
average  satisfy  the  performance  objectives  but  be  exposed  to  high- 
risk  scenarios  with  non-negligible  probabilities  [1,7,16,24,28,36], 
The  original  contributions  of  this  work  reside  in:  addressing  the 
optimal  renewable  DG  technology  selection,  sizing  and  allocation 
problem  within  a  simulation  and  multi-objective  optimization 
(MOO)  framework  that  allows  for  assessing  and  controlling  risk; 
introducing  the  conditional  value-at-risk  (CVaR)  as  a  measure  of 
the  risk  associated  to  each  objective  function  of  the  optimization 
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[37,38].  The  main  sources  of  uncertainty  are  taken  into  account 
through  the  implementation  of  a  MCS  and  OPF  (MCS-OPF)  resolu¬ 
tion  engine  nested  in  a  MOO  based  on  NSGA-II  [39],  The  aim  of 
the  MOO  is,  specifically,  the  simultaneous  minimization  of  the 
expected  global  cost  (£CS)  and  expected  energy  not  supplied 
(££NS),  and  corresponding  CVaR  values.  A  weighting  factor  /?  is 
introduced  to  leverage  the  impact  of  the  CVaR  in  the  search  of  the 
final  Pareto  optimal  renewable  DG  integration  solutions.  The 
proposed  framework  provides  a  new  spectrum  of  information 
for  well-supported  decision  making,  enabling  the  trade-off 
between  optimal  expected  performance  and  the  associated  risk 
to  achieve  it. 
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2.  Distributed  generation  network  simulation  model 

This  section  introduces  the  MCS-OPF  model,  including  the 
definition  of  the  DG  structure  and  configuration,  the  presentation 
of  the  uncertainty  sources  and  their  treatment,  the  MCS  for 
scenarios  generation  and  the  OPF  formulation  for  evaluating  the 
performance  of  the  distribution  network,  in  terms  of  the  objective 
functions  of  the  MOO  problem.  The  outputs  of  the  MCS-OPF  model 
are  the  probability  density  functions  of  the  energy  not  supplied 
( ENS )  and  the  global  cost  (Cg)  of  the  network,  and  their  respective 
CVaR  values. 

2.1.  Distributed  generation  network  structure  and  configuration 

Four  main  classes  of  components  are  considered  in  the  dis¬ 
tribution  network:  nodes,  feeders,  renewable  DG  units  and  main 
power  supply  spots  (MS).  The  nodes  can  be  understood  as  fixed 
spatial  locations  at  which  generation  units  and  loads  can  be 
allocated.  Feeders  connect  different  nodes  and  through  them  the 
power  is  distributed.  Renewable  DG  units  and  main  power  supply 
spots  are  power  sources;  in  the  case  of  electric  vehicles  and 
storage  devices  they  can  also  act  as  loads  when  they  are  in 
charging  state.  The  locations  of  the  main  supply  spots  are  fixed. 
The  MOO  aims  at  optimally  allocating  renewable  DG  units  at  the 
different  nodes.  Fig.  1  shows  an  example  of  configuration  of  a 
distribution  network  adapted  from  the  IEEE  13  nodes  test  feeder 
[40],  for  which  the  regulator,  capacitor,  switch  and  the  feeders 
with  length  equals  to  zero  are  neglected. 

Each  component  in  the  distribution  network  has  its  own 
features  and  operating  states  that  determine  its  performance. 
Assuming  stationary  conditions  of  the  operating  variables,  the 
network  operation  is  characterized  by  the  location  and  magnitude 
of  power  available,  the  loads  and  the  mechanical  states  of  the 
components,  because  degradation  or  failures  can  have  a  direct 
impact  on  the  power  availability  (in  the  DG  units,  feeders  and/or 
main  supply). 

The  renewable  DG  technologies  considered  in  this  work  include 
solar  photovoltaic  (PV),  wind  turbines  (W),  electric  vehicles  (EV) 
and  storage  devices  (batteries)  (ST).  The  power  output  of  each  of 
these  technologies  is  inherently  uncertain.  PV  and  W  generation 
are  subject  to  variability  through  their  dependence  on  environ¬ 
mental  conditions,  i.e.,  solar  irradiance  and  wind  speed.  Dis/ 
connection  and  dis/charging  patterns  in  EV  and  ST,  respectively, 
further  influence  the  uncertainty  in  the  power  outputs  from  the 
DG  units.  Also  generation  and  distribution  interruptions  caused  by 
failures  are  regarded  as  significant. 

The  following  notation  is  used  for  sets  and  subsets  of  compo¬ 
nents  in  the  distribution  networlcN— set  of  all  nodes;  MS— set  of  all 
types  of  main  supply  power  sources;  DG— set  of  all  DG  technolo¬ 
gies;  PV— set  of  all  photovoltaic  technologies;  W—  set  of  all  wind 
technologies;  £V— set  of  all  electric  vehicle  technologies;  ST—  set  of 
all  storage  technologies;  FD— set  of  all  feeders. 


Fig.  1.  Example  of  distribution  network  configuration. 


The  configurations  of  power  sources  allocated  in  the  network, 
indicating  the  size  of  power  capacity  and  the  location,  is  given  in 
matrix  form: 
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where  E—  configuration  matrix  of  type,  size  and  location  of  the 
power  sources  allocated  in  the  distribution  network;  EMS— size 
and  location  of  main  supply,  fixed  part  of  the  configuration  matrix; 
Edc— type  size  and  location  of  DG  units,  decision  variable  part  of 
the  configuration  matrix;  n— number  of  nodes  in  the  network,  INI; 
m— number  of  main  supply  type  (transformers),  IMSl;  d— number  of 
DG  technologies,  IDGI. 


fij- 


£  number  of  units  of  the  MS  type  or  DG  technology  j  allocated  at  node  1 
0  otherwise 


VieNJeMS  U  DG.f  eZ* 


(2) 


Feeders  deployment  is  described  by  the  set  of  pairs  of  nodes 
connected: 

FD=  {(1,2) . (i,i)}  V(i,ijeN  x  N,(i,i')  is  a  feeder  (3) 


Any  configuration  { E,FD }  of  power  sources  E=[Sms\Sdg]  and 
feeders  FD  of  the  distribution  network  are  affected  by  uncertainty, 
so  that  the  operation  and  performance  of  the  distribution  network 
is  strongly  dependent  on  the  network  configuration  and  scenarios. 
Furthermore,  if  the  distribution  network  acts  as  a  ‘price  taker’,  the 
variability  of  the  economic  conditions,  particularly  the  price  of  the 
energy,  is  also  an  influencing  factor  [13,19,20],  For  these  reasons,  it 
is  imperative  to  represent  and  account  for  the  uncertainties  in  the 
optimal  allocation  results  for  informed  and  conscious  decision¬ 
making. 
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2.2.  Uncertainty  modeling 
2.2.1.  Photovoltaic  generation 

PV  technology  converts  the  solar  irradiance  into  electrical 
power  through  a  set  of  solar  cells  configured  as  panels.  Commonly, 
solar  irradiance  has  been  modeled  using  probabilistic  distribu¬ 
tions,  derived  from  the  weather  historical  data  of  a  particular 
geographical  area.  The  Beta  distribution  function  [41,42]  is  used  in 
this  paper: 

,  Vs  e  [0,  l],a  >  0,3  >  0 

fpV(s)  =  {  na)m  H  (4) 

0  otherwise 


randomness  of  the  wind  speed  in  various  conditions  [1,42]: 

fw(ws)  =  2^e  - (ws/ff)2  (12) 

where  ws— wind  speed  [m/s];  fw— Rayleigh  probability  density 
function;  er— scale  parameter  of  the  Rayleigh  distribution  function. 

Then,  for  a  given  wind  speed  value,  the  power  output  of  one 
wind  turbine  can  be  determined  as  [1,41,42]: 

if  wsci  <  ws  <  wsa 

if  ws„  <  ws  <  wsco  (13) 

otherwise 


Pw(ws)  = 


1  RTD 

0 


where  s— solar  irradiance;  fpv— beta  probability  density  function;  a, 
/?— parameters  of  the  beta  probability  density  function. 

The  parameters  of  the  Beta  probability  density  function  can  be 
inferred  from  the  estimated  mean  /,/  and  standard  deviation  a  of 
the  random  variable  s  as  follows  [1]: 

,5, 


a  = 


1  —  p 


(6) 


Besides  dependence  on  solar  irradiation,  PV  depends  also  on 
the  features  of  the  solar  cells  that  constitute  the  panels  and  on 
ambient  temperature  on  site.  The  power  outputs  from  a  single 
solar  cell  is  obtained  from  the  following  equations  [41,42]: 


Tc  —  Ta +s 


fNoT-20\ 

V  0.8  ) 


(7) 


I  —  s(ISc -E  — 25))  (8) 

V  =  Voc+kvTc  (9) 


FF  = 


VmppImpp 

17  ochc 


(10) 


where  ws^— cut-in  wind  speed  [m/s];  wsa— rated  wind  speed  [m/s]; 
wsco— cut-out  wind  speed  [m/s];PPJ-D_rated  power  [kW];  Pw(ws)— 
wind  power  output  [kW], 

2.2.3.  Electric  vehicles 

In  this  work,  EV  are  considered  as  battery  electric  vehicles  with 
three  possible  operating  states:  charging,  discharging  (i.e.,  inject¬ 
ing  power  into  the  distribution  network)  and  disconnected  [43].  To 
model  their  pattern  of  operation,  they  are  considered  as  a  ‘block 
group’,  aggregating  their  single  operating  states  into  an  overall 
performance.  The  main  reasons  for  this  aggregation  are  the 
observed  nearly  stable  daily  usage  schedule  of  EV  and  the  need 
of  avoiding  the  combinatorial  explosion  of  the  model  [42], 

The  power  output  of  one  block  of  EV  is  formulated  by  assigning 
residence  time  intervals  to  each  possible  operating  state  and 
associating  them  with  the  percentage  of  trips  that  the  vehicles 
perform  by  hour  of  a  day  [43],  This  allows  approximating  the 
hourly  probability  distribution  of  the  operating  states  per  day,  as 
shown  Fig.  2.  In  a  given  (random)  scenario  of  operational  condi¬ 
tions,  the  determination  of  the  operating  state  of  a  block  of  EV,  of  a 
specific  hour  of  the  day,  is  sampled  randomly  from  the  corre¬ 
sponding  probability  distribution.  Accordingly,  the  power  output 
for  a  unit  or  block  group  of  EV  is  calculated  using  the  expressions 
below: 


(  Pdch(td)  if  op  =  discharging 

/ev(td,op)=<  Pch(^d)  if  op  =  charging  Vop  e  OPs=  {charging,  discharging,  disconnected]  (14) 

I  Pdtd(td)  if  op  =  disconnected 

!PeRYD  if  op  =  discharging 

-Prtd  if  op  =  charging  Vt  e  [0,  tRop],  op  e  OPs  =  {charging,  discharging,  disconnected)  (15) 

0  if  op  =  disconnected 


Ppv(s)  =  nceIISFF  x  V  x  I  (11) 

where  Ta— ambient  temperature  [°C];  Nor- nominal  cell  operating 
temperature  [°C];  Tc— cell  temperature  [CC];  Isc— short  circuit 
current  [A];  k,— current  temperature  coefficient  [mA/°C];  Voc— open 
circuit  voltage  [V];  kv— voltage  temperature  coefficient  [mV/°C]; 
VMPP— voltage  at  maximum  power  [V];  IMpp—  current  at  maximum 
power  [A];  FF — fill  factor;  nceiis — number  of  photovoltaic  cells; 
pPv(s)-PV  power  output  [W]. 


2.2.2.  Wind  generation 

Wind  generation  is  obtained  from  turbine-alternator  devices 
that  transform  the  kinetic  energy  of  the  wind  into  electrical  power. 
The  stochastic  behavior  of  the  wind  speed  is  commonly  repre¬ 
sented  through  probability  distribution  functions.  In  particular, 
the  Rayleigh  distribution  has  been  found  suitable  to  model  the 


where  td—  hour  of  the  day  [h] ;  tRop— residence  time  interval  for 
operating  state  op  [h] ;  fev— operating  state  probability  density 
function;  Pppd— rated  power  [kW], 


2.2.4.  Storage  devices 

Analogously  to  the  EV  case,  storage  devices  are  treated  as 
batteries.  In  reality,  these  present  two  main  operating  states, 
charging  and  discharging  [44],  However,  for  this  study  the  level 
of  charge  in  the  batteries  is  randomized  and  the  state  of  dischar¬ 
ging  is  the  only  one  that  is  allowed.  This  is  done  to  simplify  the 
behavior  of  the  batteries,  making  it  independent  on  the  previous 
state  of  charge.  The  discharging  time  interval  is  assigned  according 
to  the  relation  between  the  batteries  rated  power,  their  energy 
density  and  the  random  level  of  charge  they  present.  For  this,  the 
discharging  action  is  carried  out  at  a  rate  equal  to  the  rated  power. 
Then,  the  power  output  per  unit  of  mass  of  active  chemical  in  the 
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battery  MT  is  estimated  as  follows: 


Fig.  3.  Daily  load  profile.  Hourly  normally  distributed  load. 


/st(Qst) 


ffik  VQi'e[0,S£xMr] 
0  otherwise 


(16) 


1  RTD 


(17) 


PsC(tR)  =  Ps^D  Vt*e[0,£H  (18) 

where  (/—level  of  charge  in  the  battery  [ kj ] ;  SE— specific  energy  of 
the  active  chemical  [kj/kg];  Mj— total  mass  of  the  active  chemical 
in  the  battery  [kg];  fst— uniform  probability  density  function; 
Prtds‘— rated  power  [1<W];  t'R— discharging  time  interval  [h]. 


2.2.5.  Main  power  supply 

The  MS  spots  in  the  distribution  network  are  the  power 
stations  connected  to  the  transmission  system.  The  distribution 
transformers  are  located  on  these  spots  and  provide  the  voltage 
level  of  the  customers.  The  stochasticity  of  the  available  main 
supplies  of  power  is  represented  following  normal  distributions 
[10,45],  truncated  by  the  maximum  capacity  of  the  transformers. 

(  -Iim5)/(7m5)  v  n ms  rfi  pms  , 

/  (S)=  <  0((P™p L  ’  capj  (}g) 

0  otherwise 

where  P™— available  main  power  supply  [kW];  p.ms— Normal 
distribution  mean;  oms— Normal  distribution  standard  deviation; 
fms— Normal  probability  density  function;  P™  —  maximum  capacity 
of  the  transformer  [kW];  (J)— standard  Normal  probability  density 
function;  <P— cumulative  distribution  function  of  <p . 


where  jnc— binary  mechanical  state  variable,  /—  failure  rate  [fail- 
ures/h],  /— repair  rate  [repairs//i,],/mc— mechanical  state  probability 
mass  function. 


2.2.7.  Demand  of  power 

Overall  demands  of  power,  as  well  as  single  load  profiles  in  the 
nodes  of  the  distribution  network,  can  be  obtained  as  daily  load 
curves  in  which  to  each  hour  corresponds  one  specific  level  of 
load,  inferred  from  historical  data  [1,14,19],  In  addition,  power 
demands  profiles  can  be  considered  uncertain  following  normal 
distributions  [34], 

Within  the  proposed  modeling  framework,  the  nodal  demands 
of  power  are  defined  by  integrating  the  two  models  mentioned 
above,  i.e„  adopting  the  general  daily  load  profile  and  considering 
the  hourly  levels  of  load  as  normally  distributed.  Fig.  3  schema¬ 
tizes  the  previous  assumption  for  a  generic  node  i. 

In  this  manner,  the  nodal  demand  of  power  is  deducted  from 
the  overall  demand  in  the  network,  and  modeled  as: 


fu(h,  td) 


'l-<P(-p,(td)/CT,(td)) 


V  i  e  N.  Li  e  [0,  oo] 


0  otherwise 


(22) 


where  td— hour  of  the  day  [h],  L— power  demand  in  node  i  [kW], 
/q— normal  distribution  mean  of  power  demand  in  node  i,  a— 
normal  distribution  standard  deviation  of  power  demand  in  node 
i,  fLi— normal  probability  density  function  of  power  demand  in 
node  i. 


2.3.  Monte  Carlo  simulation 


2.2.6.  Mechanical  states  of  the  components 

Renewable  DG  units,  MS  spots  and  feeders  are  subject  to  wearing 
and  degradation  processes.  These  processes  can  trigger  unexpected 
events,  even  failures,  interrupting  or  reducing  the  specific  function¬ 
ality  of  each  component.  Frequently,  the  stochastic  behavior  of 
failures,  repairs  and  maintenance  actions  is  modeled  using  Markov 
models  [28,42],  In  this  work,  a  two-state  model  is  implemented  in 
which  the  components  can  be  in  the  mutually  exclusive  states: 
available  to  operate  and  under  repair  (failure  state).  Assuming  the 
duration  of  each  state  as  exponentially  distributed,  the  mechanical 
state  of  a  component  can  be  randomly  generated  as  follows: 

ft  if  the  component  is  available  to  operate  „ 

mc={  Vcomponente  I^.FD) 

\  0  otherwise 

(20) 


fmcimc)  = 


(1  -me)/  +  mc/ 

FTF 


Vmc  e  [0, 1) 


(21) 


Most  of  the  techniques  used  for  evaluating  the  performance  of 
renewable  DG-integrated  distribution  networks  are  of  two  classes: 
analytical  methods  and  MCS  [28].  The  implementation  of  analy¬ 
tical  methods  is  always  preferable,  in  theory,  because  of  the 
possibility  of  achieving  closed  exact  solutions,  but  in  practice;  it 
often  requires  strongly  simplifying  assumptions  that  may  lead  to 
unrealistic  results:  power  network  applications  exist  but  for  non- 
fluctuating  or  non-intermittent  generation  and/or  load  profiles, 
and  low  dimensionality  of  the  network,  gaining  traceability  with 
reduced  computational  efforts  [32],  Different,  MCS  techniques 
allow  considering  more  realistic  models  that  analytical  methods 
do,  because  simplifying  assumptions  are  not  necessary  to  solve  the 
model,  since  de  facto  the  model  is  not  solved  but  simulated  and  the 
quantities  of  interest  are  estimated  from  the  statistics  of  the  virtual 
simulation  runs  [46],  For  this  reason  MCS  is  quite  adequate  for 
application  on  the  analysis  of  distribution  networks  with  significant 
randomness  or  variability  in  the  sources  of  power  supply  and 
loads,  failure  occurrence  and  strong  dependence  on  the  power 
flows  as  a  consequence  of  congestion  conditions  in  the  feeders,  etc. 
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Fig.  4.  Example  of  the  matrix  form  construction  of  a  DG-integrated  network 
(A)  and  schema  of  the  operating  state  definition  from  the  sampled  variables  (B). 


[3,31,33,41,42,47];  the  price  to  pay  for  this  is  the  possibly  consider¬ 
able  increment  in  the  use  of  computational  resources,  and  various 
methods  exist  to  tackle  this  problem  [46], 

Given  the  multiple  sources  of  uncertainties  considered  in  the 
proposed  framework  and  the  proven  advantages  of  MCS  for 
adequacy  assessment  of  power  distribution  networks  with  uncer¬ 
tainties  [3,31,33,41,42,47],  we  adopt  a  non-sequential  MCS  to 
emulate  the  operation  of  a  distribution  network,  sampling  the 
uncertain  variables  without  considering  their  time  dependence,  so 
as  to  reduce  the  computational  problem. 

For  a  given  structure  and  configuration  of  the  distribution 
network  { E,FD },  i.e„  for  the  fixed  EMS  and  FD  deployments  and 
the  proposed  renewable  DG  integration  plan  denoted  by  EDG,  each 
uncertain  variable  is  randomly  sampled.  The  set  9  of  sampled 
variables  constitutes  an  operational  scenario,  in  correspondence  of 
which  the  distribution  network  operation  is  modeled  by  OPF  and 
its  performance  evaluated.  The  two  inputs  to  the  OPF  model  are 
the  network  configuration  { E,FD }  and  the  operational  conditions 
scenario  9. 

9  =  tvs,-,  Qfj,  mcjj,  mc^j]  Vi,i'  e  N.j  e  MS  u  DG,  (i,  i')  e  FD 

(23) 


the  generation  unit  j,  formulated  in  Eqs.  (24)  and  (25). 

PL,  =  kimcuGj<9) 


(24) 


Gj(9)  =  < 


ptns;  8 

if  j  e  MS 

P^V(sf) 

if  j  e  PV 

Ppwsf) 

if  j  e  W 

Pfv(°P?j ) 

if  j  e  EV 

Pf(Qf) 

if  j  s  ST 

(25) 


In  the  proposed  non-sequential  MCS  procedure,  the  intermit- 
tency  in  the  solar  irradiation  is  taken  into  account  defining  a  night 
interval  between  22.00  and  06.00  h,  i.e„  if  the  value  of  the  hour  of 
the  day  td  (h),  sampled  from  a  discrete  uniform  distribution  U 
(1,24),  falls  in  the  night  interval,  there  is  no  solar  irradiation. 
Regarding  the  wind  speed,  its  variability  is  considered  by  sampling 
positive  values  from  a  Rayleigh  probability  density  function  fitted 
on  historical  data  and  whose  parameters  as  such  that  the  prob¬ 
ability  of  absence  of  wind  is  zero.  Since  it  is  not  reasonable  to  force 
the  historical  profile  of  the  wind  speed  to  follow  a  distribution  that 
admits  intermittency,  a  common  alternative  technique  is  to  model 
the  wind  by  a  Markov  Chain.  Indeed,  it  is  possible  to  accurately 
represent  the  wind  speed  by  a  stationary  Markov  process  if  the 
historical  profile  of  wind  speed  data  is  sufficiently  large  e.g.,  years 
[28],  The  intermittency  is,  then,  represented  by  the  first  state  of 
the  chain  with  wind  speed  equals  to  zero,  and  the  sampling  of  the 
wind  speed  states  in  the  non-sequential  MCS  of  the  proposed 
framework,  can  be  performed  using  the  steady-state  probabilities 
of  the  Markov  Chain. 

An  important  issue  in  modeling  the  operation  of  power 
systems  is  how  to  represent  the  evolution  of  uncertain  operating 
conditions,  such  as  solar  irradiation,  wind  speed,  load  profiles, 
energy  prices,  among  others.  As  an  example,  the  load  forecast 
implies  the  prediction  of  future  power  demands  given  specific 
previous  conditions.  Therefore,  to  consider  load  forecast  uncer¬ 
tainty  within  the  proposed  MCS  framework,  it  would  be  necessary 
to  change  to  a  sequential  simulation  model,  in  which  the  uncertain 
renewable  energy  resources,  main  power  supply  and  loads  must 
be  sampled  at  each  time  step.  In  particular,  load  forecast  uncer¬ 
tainty  can  be  integrated  properly  building  consecutive  load  sce¬ 
narios  and  assigning  corresponding  probabilities  of  occurrence  as 
presented  by  [7,48],  Another  interesting  approach  for  load  forecast 
uncertainty  modelling  is  the  geometric  Brownian  motion  (CBM) 
stochastic  process  [31,49], 

2.4.  Optimal  power  flow 


where  td— hour  of  the  day  [h],  randomly  sampled  from  a  discrete 
uniform  distribution  U(l,24). 

Fig.  4  shows  an  example  of  the  matrix  form  construction  of  the 
DG-integrated  distribution  network,  considering  a  simple  case  of 
n= 3  nodes.  The  network  contains  one  MS  spot  at  node  i= 1, 
defining  the  fixed  part  EMS  of  the  configuration  matrix,  whereas, 
the  decision  variable  EDG  proposes  a  renewable  DG  integration 
plan  EDC  that  built  from  the  number  of  units  £  of  each  DG 
technology  allocated.  In  this  way,  the  network  configuration 
{.E.FD}  is  composed  by  the  matrix  E=[EMS  I  EDG]  and  the 
deployment  of  feeders.  Then,  given  the  spatial  representation  {S', 
FD },  the  sampling  of  the  scenario  9  determines  the  operational 
conditions  to  perform  power  flow  analysis,  i.e.,  distribute  the 
power  available  P®0  to  supply  appropriately  the  demands  Lt. 

The  available  power  in  the  power  source  type  j  at  node  i,  Pf;a .  is 
function  of  the  number  of  units  allocated  £y,  the  mechanical  state 
mCy  and  the  specific  unitary  power  output  function  G/associated  to 


Power  flow  analysis  is  performed  by  DC  OPF  [50]  which  takes 
into  account  the  active  power  flows,  neglecting  power  losses,  and 
assumes  a  constant  value  of  the  voltage  throughout  the  network. 
This  allows  transforming  to  linear  the  classic  non-linear  power 
flow  formulation,  gaining  simplicity  and  computational  tractabil- 
ity.  For  this  reason,  DC  power  flow  is  often  used  in  techno- 
economic  analysis  of  power  systems,  more  frequently  in  transmis¬ 
sion  [50,51]  but  also  in  distribution  networks  [51]. 

The  DC  power  flow  generic  formulation  is: 

pi=  2  Vi,t'eJV,(i,i')eFD  (26) 

i'  eN  '  ' 

I,  (Pa  -  Li- pi)  =  0  WieN  (27) 

where,  P— active  power  leaving  node  i  [1<W[;  Bir  susceptance  of 
the  feeder  (i,f)  [1/Q[;  <5j— voltage  angle  at  node  i;  PCl—  active  power 
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injected  or  generated  at  node  i  [1<W];  L,— load  at  node  i  [kWj.The 
assumptions  are: 

•  the  difference  between  voltage  angles  are  small,  i.e.,  sin(AS) 
«  8,  cos(A8)  1 

•  the  feeders  resistance  are  neglected,  i.e.,  R«X,  which  implies 
that  power  losses  in  the  feeder  are  also  neglected 

•  the  voltage  profile  is  flat  (constant  V,  set  to  1  p.u.) 


Then,  for  a  given  configuration  {S.FD}  and  operational  scenario 
$  the  formulation  of  the  OPF  problem  is: 


min  Cn£»(P9u)  =  £  2  C0&MJP£  ts 

ieNj  eMSuDG  1 


(28) 


s.t. 


Li-  X  X™c?ji,B(i/)(S*-S?)-LS?  =  0  Vi,i'EN,(i,i')eFD 

j  e  MSUDG  f  eN 


Pcu^PL,  Vi  g  N,j  e  MS  U  DG 


(29) 

(30) 


0  <  P9U.  V  i  e  N,j  e  MS  U  DG 


(31) 


2.5.2.  Global  cost 

The  Cg  of  the  distribution  network  is  formed  by  two  terms, 
fixed  and  variable  costs.  The  former  term  includes  those  costs  paid 
at  the  beginning  of  the  operation  after  the  installation  of  the  DG 
(conception  of  SDG).  They  are  the  investment-installation  cost  and 
the  operation-maintenance  fixed  cost.  The  variable  term  refers  to 
the  operating  and  maintenance  costs.  Note  that  these  costs  are 
dependent  on  the  power  generation  and  supply,  which  are  a  direct 
output  of  the  OPF  (Eq.  (28)).  In  addition,  this  term  considers 
revenues  associated  to  the  renewable  sources  incentives.  Consid¬ 
ering  the  distribution  network  as  a  ‘price  taker’  entity,  the  profits 
depend  on  the  value  of  the  energy  price  that  is  correlated  with  the 
total  load  in  the  network.  Three  different  ranges  of  load  are 
considered  for  the  daily  profile.  For  each  range,  a  correlation  value 
of  energy  price  is  considered  as  shown  in  Fig.  5(A). 

In  Fig.  5(B)  the  correlation  between  energy  price  and  total  load 
is  presented  as  the  proportion  of  their  maximum  values.  As  an 
intermediate  approximation  of  existing  studies  (e.g.,  [13,19,20]), 
the  line  with  square-markers  represents  the  proportional  correla¬ 
tion  used  in  this  study,  which  can  be  expressed  as: 

ep  =  ep(l^--0.38(i^))2  +  1.38igf)j  (36) 


mc(if)B(i,?)(fi?  x  Amp(i/)  V  i,  i  e  N,  ( i ,  i')  e  FD  (32) 

-  -S9)<Vx  Ampan  V  I,  i'  e  N,  (i,  i')  e  FD  (33) 

where,  ^—duration  of  the  scenario  [h];  Cg^— operating  and 
maintenance  costs  of  the  total  power  supply  and  generation 
[ $ ] : Go&m"  — operating  and  maintenance  variable  costs  of  the  power 
source  j  [$/l<Wh[;  mrfj,.—  mechanical  state  of  the  feeder  (i,i'); 
Btf— susceptance  of  the  feeder  (i,i'),  [1/Q];  rncfj,  mechanical  state 
of  the  power  source  j  at  node  i;  P9a.—  available  power  in  the  source 
j  at  node  i  [kW[;  P9U power  produced  by  source)  at  node  i  [kW[; 
LSf— load  shedding  at  node  i  [1<W[;  V— nominal  voltage  of  the 
network  [kV[;  Amp^— ampacity  of  the  feeder  (i,i'),  [A]. 

The  load  shedding  in  the  node  i,  LSif  is  defined  as  the  amount  of 
load(s)  disconnected  in  node  i  to  alleviate  overloaded  feeders  and/ 
or  balance  the  demand  of  power  with  the  available  power  supply 
[52], 

The  OPF  objective  is  the  minimization  of  the  operating  and 
maintenance  costs  associated  to  the  generation  of  power  for  a 
given  scenario  8  of  duration  Is.  Eq.  (29)  corresponds  to  the  power 
balance  equation  at  node  i,  while  Eqs.  (30)  and  (31 )  are  the  bounds 
of  the  power  generation  and  Eqs.  (32)  and  (33)  account  for  the 
technical  limits  of  the  feeders. 

2.5.  Performance  indicators 

Given  a  set  T  of  ns  sampled  operational  scenarios  8?,  €  e  {1,. .  ,,ns}, 
the  OPF  is  solved  for  each  scenario  8e  e  T,  giving  in  output  the  values 
of  ENS  and  global  cost. 

2.5.3.  Energy  not  supplied 

ENS  is  a  common  index  for  reliability  evaluation  in  power 
systems  [1,10,11,48,49,52-55].  In  the  present  work,  its  value  is 
obtained  directly  from  the  OPF  output  in  the  form  of  the  aggrega¬ 
tion  of  all-nodal  load  sheddings  per  scenario  8?: 

ENS9'  =  2  FS9'  x  ts  V.9,  e  r 

i  e  N 

ENSr  =  {ENSs\..„  ENS9' ENS9"1 } 


(34) 

(35) 


Thereby,  the  global  cost  function  for  a  scenario  i9^  is  given  by: 


qr=  2  z  (Cfnvj  ■ 

i  e  Nj  e  DG 


-C, 


-(inc+ep(L9'))  2  2 


(37) 


i  e  Nj  e  DG 


=  {C9\...,Cg',--,Cgns}  (38) 

where  CinVj— investment  cost  of  the  DG  technology  j  [$];  C0&Mf  — 
operating  and  maintenance  fixed  costs  of  the  DG  technology 
j  [$];  t!1— horizon  of  analysis  [h[;  inc— incentive  for  generation  from 
renewable  sources  [$/kWh[;  ep— energy  price  [$/l<Wh[;  C9'— 
global  cost  [$]. 


2.5.3.  Risk 

In  [38],  the  importance  of  measuring  risk  when  optimizing 
under  uncertainty  and  including  it  as  part  of  the  objective  function 
(s)  or  constraints  is  emphasized.  The  proposed  MOO  framework 
introduces  the  CVaR  as  a  coherent  measure  of  the  risk  associated 
to  the  objective  functions  of  interest.  The  CVaR  has  been  broadly 
used  in  financial  portfolio  optimization  either  to  reduce  or  mini¬ 
mize  the  probability  of  incurring  in  large  losses  [37,38].  This  risk 
measurement  allows  evaluating  how  ‘risky’  is  the  selection  of  a 
solution  leading  to  a  determined  value  of  expected  losses. 

We  can  consider  a  fixed  configuration  of  the  distribution 
network  {S.FD}  including  the  integration  of  DG  units  as  a 
‘portfolio’.  The  assessed  expectations  of  ENS}  and  cj,  found  from 
the  MCS-OPF  applied  to  the  set  of  scenarios  T,  are  estimations  of 
the  ‘losses’;  then,  CVaR(ENSr )  and  CVaR(Cg)  represent  the  risk 
associated  to  the  solutions  with  these  expectations. 

The  definition  of  CVaR  for  continuous  and  discrete  general  loss 
functions  is  given  in  detail  in  [38],  Here  a  simplified  and  intuitive 
manner  to  understand  the  CVaR  definition  and  its  derivation 
according  to  [56]  is  presented. 

As  shown  in  Fig.  6(A),  for  a  discrete  approximation  of  the 
probability  of  the  losses,  given  a  confidence  level  or  a-percentile, 
the  value-at-risk  VaRa  represents  the  smallest  value  of  losses  for 
which  the  probability  that  the  losses  do  not  exceed  the  value  of 
VaRa  is  greater  than  or  equal  to  a.  Thus,  from  the  cumulative 
distribution  function  F(losses)  is  possible  to  construct  the  a-tail 
cumulative  distribution  function  Fa(losses )  for  the  losses,  such  that 
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Lj{td)/LTh 


Fig.  5.  Example  of  load  ranges  definition  for  a  generic  daily  load  profile  (A)  and  correlation  energy  price-total  load  (B)  [13,19,20]. 
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Fig.  6.  Graphic  representation  of  the  CVaR. 


(Fig.  6(B)): 


3.1.  MOO  problem  formulation 


Fa{losses) 


F(/osses^  a  jf  Yapa  <  losses 

0  otherwise 


(39) 


The  a-tail  cumulative  distribution  function  represents  the  risk 
‘beyond  the  VaR’  and  its  mean  value  corresponds  to  the  CVaRa. 

Among  other  risk  measures,  the  CVaR  has  been  commonly  used 
to  assess  the  financial  impact  associated  to  different  sources  of 
uncertainty  on  electricity  markets  behavior.  Some  interesting 
approaches  in  the  use  of  diverse  risk  measures  for  electricity 
markets  modelling  can  be  found  in  [49,57,58] 


Considering  a  set  of  randomly  generated  scenarios  Y,  the 
optimization  problem  is  formulated  as  follows: 

min  fA=PECrg+(\-P)CVaRa(Crg)  (40) 


min  f2  =  f)EENSr  +  (1  -  fi)CV  aRa(ENSr ) 


(41) 


s.t. 


£  number  of  units  of  the  MS  type  or  DG  technology  j  allocated  at  node  / 
0  otherwise 


Vi  e  N,j  g  MS  U  DG,  £  e  Z* 


(2) 


3.  DC  units  selection,  sizing  and  allocation 


z  z 

i  e  Nj  e  DG 


lij(Cj„v.  +  C0&Mr )  <  BCT 


(42) 


This  section  presents  the  general  formulation  of  the  MOO 
problem  considered  previously.  As  introduced,  the  practical  aim 
of  the  MOO  is  to  find  the  optimal  integration  of  DG  in  terms  of 
selection,  sizing  and  allocation  of  the  different  renewable  genera¬ 
tion  units  (including  EV  and  ST).  The  corresponding  decision 
variables  are  contained  in  SDC  of  the  configuration  matrix  3. 

The  MOO  problem  consists  in  the  concurrent  minimization  of 
the  two  objective  functions  measuring  the  Cg  and  ENS,  and  their 
associated  risk.  Specifically,  their  expected  values  and  their  CVaR 
values  are  combined,  weighted  by  a  factor  [0,1],  which  allows 
modulating  the  expected  performance  of  the  distribution  network 
and  its  associated  risk. 


Z ‘Jij  —  TJ  VjeDG  / 43> 

ieN  ’ 

OPF({S,  FD),T)  (28)— (33) 

where  ECg  and  EENS  denote  the  expected  values  of  Cg  and  ENS, 
respectively. 

The  meaning  of  each  constraint  is,  (2)— the  decision  variable 
is  a  non-negative  integer  number;  (42)— the  total  costs  of  invest¬ 
ment  and  fixed  operation  and  maintenance  of  the  DG  units  must 
be  less  or  equal  to  the  available  budget  BCT,  (43)— the  total 
number  of  DG  units  to  allocate  of  each  technology  j  must  be  less 
or  equal  to  the  maximum  number  of  units  available  Tj  to  be 
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NSGA-II 


Fig.  7.  Flow  chart  of  NSGA-II  MCS-OPF  MOO  framework. 
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integrated;  (28)-(33)— all  the  equations  of  OPF  must  be  satisfied 
for  all  scenarios  in  Y. 

Constraint  (43)  can  be  translated  into  maximum  allowed 
penetration  factor  PF°  jX.  °f  each  DG  technology  j.  Defining  PF  as 
‘the  output  active  power  of  total  capacity  of  DG  divided  by  the 
total  network  load’  [59],  constraint  (43)  can  be  rewritten  as 
follows: 


I  €ijEPj 

i'eN 

ELt 


TjEPj 

;Jar  vjeDG 


(44) 


Fig.  8.  Radial  11 -nodes  distribution  network. 


Table  1 

Feeders  characteristic  and  technical  data  [40], 


Type 

Node  i 

Node  i' 

Length  [km] 

X  [n/km] 

Am p  [A] 

12 

1 

2 

0.61 

0.37 

365 

72 

2 

3 

0.15 

0.47 

170 

73 

2 

4 

0.15 

0.56 

115 

72 

2 

6 

0.61 

0.37 

365 

73 

4 

5 

0.09 

0.56 

115 

76 

6 

7 

0.15 

0.25 

165 

74 

6 

8 

0.09 

0.56 

115 

72 

6 

11 

0.31 

0.37 

365 

75 

8 

9 

0.09 

0.56 

115 

77 

8 

10 

0.24 

0.32 

115 

Table  2 

Main  power  supply  parameters. 

Node  i  P^p  [kW] 

Normal  distribution  parameters 

V™  ams 

1  1600 

1200  27.5 

where  £  £y— is  the  total  number  of  units  of  DG  technology  j 

i'eN 

integrated  in  the  network;  £P°C— is  the  expected  power  output  of 
one  unit  of  DG  technology  j  [kW[;  ELT— is  the  expected  total 
load  [kW[. 

The  MOO  optimization  problem  is  non-linear  and  non-convex, 
i.e„  a  non-convex  mixed-integer  non-linear  problem  or  non- 
convex  MINLP.  It  is  non-linear  because  the  objective  functions 
given  by  Eqs.  (40)  and  (41 )  cannot  be  written  in  the  canonical  form 
of  a  linear  program,  i.e„  CTX,  where  C  a  vector  of  known 
coefficients  and  X  the  decision  vector.  In  the  present  case,  the 
decision  matrix  EDC  enters  the  MCS-OPF  flow  simulation  to  obtain 
the  probability  mass  functions  of  Cg  and  ENS  and,  then,  the 
objective  functions  are  formed  from  the  corresponding  expected 
and  CVaR  values.  Thus,  the  operations  applied  on  EDC  through 
MCS-OPF,  expectation  and  CVaR  cannot  not  be  represented  as  the 
product  CTEDC.  The  problem  is  non-convex  because  the  decision 
matrices  EDC  are  integer-valued  (constraint  (2))  and,  as  it  is 
known,  the  set  of  non-negative  integers  is  non-convex. 

Given  the  class  of  optimization  problems  in  the  proposed 
framework  (non-convex  MINLP),  it  is  most  likely  to  have  multiple 
local  minima.  Moreover,  the  dimension  of  the  distribution  net¬ 
work  can  lead  to  a  combinatorial  explosion  of  the  feasible  space  of 
the  decision  matrices  EDC  [7,10],  incrementing  the  number  of 
possible  local  minima  and  hindering  the  possibility  of  benchmark¬ 
ing  the  optimal  solutions  obtained.  However,  an  approximated  but 
straightforward  alternative  is  to  perform  several  realizations  of  the 
framework  obtaining  different  optimal  solutions  under  the  same 
optimization  and  simulation  conditions  (parameters)  and,  thus, 
compare  them  regarding  the  optimal  decision  matrices  and  their 
associated  value  of  the  objective  functions. 


Table  3 

Parameters  of  PV,  W,  EV  and  ST  technologies  [11,13,42]. 


PV 

W 

Beta  distribution  a 

0.26 

Rayleigh  distribution  o 

7.96 

Beta  distribution  p 

0.73 

Prtd  [kW] 

50 

Ta  1  C] 

30 

wsci  [m/s] 

3.8 

Not  l  C] 

43 

wsa  [m/s] 

9.5 

4c  [A] 

1.8 

wsco  [m/s] 

23.8 

ki  [mA/’C] 

1.4 

EV 

Voc  IV] 

55.5 

PZtd  [kW] 

6.3 

k„  [mV/°C] 

194 

ST 

Vmpp  [V] 

38 

PStd  [kW] 

0.275 

Impp  [A] 

1.32 

SE  [kj/kg] 

0.042 

A 


B 


Fig.  9.  Mean  (A)  and  variance  (B)  values  of  nodal  power  demand  daily  profiles. 
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This  process  was  performed  for  the  proposed  case  study. 
Indeed,  the  optimal  decision  matrices  EDC  are  different  in  all  the 
cases,  when  the  optimization  and  simulation  framework  is  per¬ 
formed  under  the  same  conditions  but,  nonetheless,  practically  the 
same  Pareto  optimal  values  of  ECg  and  EENS  are  eventually 
obtained.  This  reflects  that  equally  expected  performances  ( ECg , 
EENS)  can  be  obtained  for  different  EDG  considering  the  large 
amount  of  feasible  combinations,  which  is  what  is  of  interest  for 
practical  applications. 

3.2.  NSGA-II  with  nested  MCS-OPF 

The  combinatorial  MOO  problem  under  uncertainties  is  solved 
by  the  NSGA-II  algorithm  [39],  in  which  the  evaluation  of  the 
objective  functions  is  performed  by  the  developed  MCS-OPF.  The 
NSGA-II  is  one  of  the  most  efficient  evolutionary  algorithms  to 
solve  MOO  problems  [60],  The  extension  to  MOO  entails  the 
integration  of  Pareto  optimality  concepts.  In  general  terms,  solving 
a  MOO  problem  of  the  form: 

min  [/j(X),/2(X),  ...Jk(X)} 
subject  toXe  A 

with  at  least  two  conflicting  objectives  functions  (/):  5Hn  — >  5H) 
implies  to  find,  within  a  set  of  acceptable  solutions  that  belong 
to  the  non-empty  feasible  region  A  <=  1R",  the  decision  vectors 
XeA  that  satisfy  the  following  [61]: 

-.3XeA//((X)  s/,(X')  Vi=  1 . k  and  /,(X)  </i(X')  for  at  least  one  i 

u- 

/(X)  </(X')  i.e./(X)  dominates /(X')  (46) 

X  is  called  a  Pareto  optimal  solution  and  the  Pareto  front  PF  is 
defined  as  {f[X)  e  5R/X  is  Pareto  optimal  solution }. 

The  process  of  searching  the  non-dominated  solutions  set  PF, 
carried  out  by  the  NSGA-II  MCS-OPF,  can  be  summarized  as  shown 
in  Fig.  7. 

The  interested  reader  can  consult  [62-64]  to  compare  the 
proposed  framework  to  alternative  MOO  analytical  approaches  in 
energy  applications. 


4.  Case  study 

We  consider  a  distribution  network  adapted  from  the  IEEE  13 
nodes  test  feeder  [40,65],  The  spatial  structure  of  the  network  has 
not  been  altered  but  we  neglect  the  regulator,  capacitor  and 
switch,  and  remove  the  feeders  of  zero  length.  The  network  is 
chosen  purposely  small,  but  with  all  relevant  characteristics  for 
the  analysis,  e.g.,  comparatively  low  and  high  spot  and  distributed 
load  values  and  the  presence  of  a  power  supply  spot  [65],  The 
original  IEEE  13  nodes  test  feeder  is  dimensioned  such  that  the  total 
power  demand  is  satisfied  without  lines  overloading.  We  modify  it 
so  that  it  becomes  of  interest  to  consider  the  integration  of 
renewable  DG  units.  Specifically,  the  location  and  values  of  some 
of  the  load  spots  and  the  ampacity  values  of  some  feeders  have 
been  modified  in  order  to  generate  conditions  of  power  congestion 
of  the  lines,  leading  to  shortages  of  power  supply  to  specific 
portions  of  the  network. 

4.1.  Distribution  network  description 

The  distribution  network  presents  a  radial  structure  of  n  =  ll 
nodes  and/d=(n-l)=10  feeders,  as  shown  in  Fig.  8.  The  nominal 
voltage  is  V=4.16  [kV],  constant  for  the  resolution  of  the  DC 
optimal  power  flow  problem. 

Table  1  contains  the  technical  characteristics  of  the  different 
types  of  feeders  considered:  specifically,  the  indexes  of  the  pairs  of 


nodes  that  are  connected  by  each  feeder  of  the  network,  their 
length,  reactance  X  and  their  ampacity  Amp. 

Concerning  the  main  power  supply  spot,  the  maximum  active 
power  capacity  of  the  transformer  and  the  parameters  of  the  normal 
distribution  that  describe  its  variability  are  given  in  Table  2. 

The  nodal  power  demands  are  reported  as  daily  profiles, 
normally  distributed  on  each  hour.  The  mean  p  and  variance  a 
values  of  the  nodal  daily  profiles  of  the  power  demands  are  shown 
in  Fig.  9(A)  and  (B),  respectively. 

The  technical  parameters  of  the  four  different  types  of  DG 
technologies  available  to  be  integrated  into  the  distribution  net¬ 
work  (PV,  W,  EV  and  ST)  are  given  in  Table  3.  The  values  of  the 
parameters  of  the  Beta  and  Rayleigh  distributions  describing  the 
variability  of  the  solar  irradiation  and  wind  speed,  are  assumed 
constant  in  the  whole  network,  i.e.,  the  region  of  distribution  is 
such  that  the  weather  conditions  are  the  same  for  all  nodes. 

The  hourly  per  day  operating  states  probability  profile  of  the  EV 
is  presented  in  Fig.  10  and  failures  and  repair  rates  of  the 
components  of  the  distribution  network  are  provided  in  Table  4. 

The  values  of  the  investment  (Cinv)  and  fixed  and  variable 
Operational  and  Maintenance  (C0&M/  and  C0&Mv)  costs  of  the  MS 
and  DG  units  are  reported  in  Table  5.  Consistently  with  the 
constraints  (42)  and  (43)  of  the  MOO  problem,  the  total  invest¬ 
ment  associated  to  a  decision  variable  EDG  (proposed  by  the 
NSGA-II)  must  be  less  than  or  equal  to  the  limit  budget;  which  is 
set  to  BGT=4500000  [$],  and  the  total  number  of  units  of  each 
type  of  DG  (following  the  order  [PV,  W,  EV,  ST])  must  be  less  than 


Fig.  10.  Hourly  per  day  probability  data  of  EV  operating  states. 


Table  4 

Failure  rates  of  feeders,  MS  and  DG  units  [11,13,42,66]. 


Type 

XF  [failures/h] 

XR  [repairs/h] 

MSuDG 

FD 

MSuDG 

FD 

MSuDG 

FD 

MS 

71 

3.33E-04 

3.33E-04 

0.021 

0.198 

PV 

72 

4.05E-04 

4.05E-04 

0.013 

0.162 

w 

73 

3.55E-04 

3.55E-04 

0.015 

0.185 

EV 

74 

3.55E-04 

3.55E-04 

0.105 

0.185 

ST 

75 

3.55E-04 

3.55E-04 

0.073 

0.185 

- 

76 

- 

4.00E-04 

- 

0.164 

- 

77 

- 

3.55E-04 
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or  equal  to  r=[  15000,  5,  200,  8000].  The  value  of  the  incentive  for 
renewable  kW  h  supplied  is  taken  as  0.024  [S/kW  h]  [34],  The 
maximum  value  of  the  energy  price  eph  is  0.11  [$/kWh]  [19,20], 
Concerning  the  calculation  of  the  CVaR,  the  alpha-percentile  is 
taken  as  a=0.80. 

Five  optimizations  runs  of  the  NSGA-II  with  the  nested  MCS-PF 
algorithm  have  been  performed,  each  one  with  a  different  value  of 
the  weight  parameter  f)  e  {1,  0.75,  0.5,  0.25,  0},  to  analyze  different 
tradeoffs  between  optimal  average  performance  and  risk.  From 
Eqs.  (40)  and  (41),  note  that  the  value  /?=1  corresponds  to 
optimizing  only  the  expected  values  of  ENS  and  Cg,  whereas  /?= 0 
corresponds  to  the  opposite  extreme  case  of  optimizing  only  the 
CVaR  values.  Each  NSGA-II  run  is  set  to  perform  g=  300  genera¬ 
tions  over  a  population  of  sz=  100  chromosomes  and,  for  the 
reproduction,  the  single-point  crossover  and  mutation  genetic 
operators  are  used.  The  crossover  probability  is  pco=  1,  whereas 
the  mutation  probability  is  pmu  =  0A ;  the  mutation  can  occur 
simultaneously  in  any  bit  of  the  chromosome. 

Finally,  sn =250  random  scenarios  are  simulated  by  the  MCS-OPF 
with  time  step  £*=1  [h].  Over  an  horizon  of  analysis  of  10  years 
(th=  87600  [h]),  in  which  the  investment  and  fixed  costs  are  prorated 
hourly. 


Table  5 

Investment,  fixed  O&M  and  variable  O&M  costs  of  MS  and  DG  [27,34,66], 


Type 

Qnv  +  ^o&A/ 

Q)&mv  I^W  h] 

MS 

_ 

1.45E-01 

PV 

48 

3.76E-05 

W 

113750 

3.90E-02 

EV 

17000 

2.20E-02 

ST 

135.15 

4.62E-05 

4.2.  Results  and  discussion 

The  Pareto  fronts  resulting  from  the  NSGA-II  MCS-OPF  are 
presented  in  Fig.  11  for  the  different  values  of  f).  The  ‘last  generation' 
population  is  shown  and  the  non-dominated  solutions  are  marked 
in  bold. 

Each  non-dominated  solution  in  the  different  Pareto  fronts 
corresponds  to  an  optimal  decision  matrix  EDC  for  the  sizing  and 
allocation  of  DG,  i.e.,  an  optimal  DG-integrated  network  config¬ 
uration  { 3, ED }  where  5=  [SMSI.3DC]. 

In  the  Pareto  fronts  obtained,  we  look  of  three  representative  non- 
dominated  solutions  for  the  analysis:  those  with  minimum  values  of 
the  objective  functions /i  and  f2  independently  and  in^2, 

respectively)  and  an  intermediate  solution  at  the  ‘elbow’  of  the  Pareto 
front.  Table  6  presents  the  values  of  the  objective  functions,  EENS,  ECg 
and  their  respective  CVaR  values  for  the  selected  solutions.  The  EENS, 
ECg  and  CVaR  values  of  the  case  in  which  no  DG  is  integrated  in  the 
network  (MS  case)  is  also  reported. 

Fig.  12  shows  a  bubble  plot  representation  of  the  selected 
optimal  solutions.  The  axes  report  the  EENS  and  ECg  values  while 
the  diameters  of  the  bubbles  are  proportional  to  their  respective 
CVaR  values.  The  MS  case  is  also  plotted. 

From  Table  6  and  Fig.  12  it  can  be  seen  that,  the  MS  case  has  an 
expected  performance  {EENS=  1109.21  [kWh]  and  £Cg=  170.27  [$]) 
inferior  (high  EENS  and  ECg )  to  any  case  for  which  DG  is  optimally 
integrated.  Furthermore,  the  CVaR(ENS)=  1656.53  [kWh]  for  the  MS 
case  is  the  highest,  indicating  the  high  risk  of  actually  achieving  the 
expected  performance  of  energy  not  supplied.  This  confirms  that  DG  is 
capable  of  providing  a  gain  of  reliability  of  power  supply  and  economic 
benefits,  the  risk  of  falling  in  scenarios  of  large  amounts  of  energy  not 
supplied  being  reduced. 

Comparing  among  the  selected  optimal  DG-integrated  networks,  in 
general  the  expected  performances  of  EENS  and  ECg  are  progressively 
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Fig.  11.  Pareto  fronts  for  different  values  of  p. 
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Table  6 

Objective  functions:  expected  and  CVaR  values  of  selected  Pareto  front  solutions. 


P 

/,  [kWh] 

f2  [kW  h] 

EENS  [kW  h] 

CVaR(ENS)  [kW  h] 

ECg  [$] 

CVaR(Cg)  [$] 

MS 

_ 

_ 

_ 

1109.21 

1656.53 

170.27 

179.24 

-DG 

“min/j 

1.00 

666.95 

160.91 

666.95 

1093.12 

160.91 

185.11 

-DG 

671.05 

150.83 

671.05 

1185.53 

150.83 

179.47 

-DG 
“min  f2 

726.57 

148.68 

726.57 

1279.37 

148.68 

178.23 

-DC 
“min  fx 

0.75 

797.07 

166.41 

677.74 

1155.11 

160.68 

183.62 

-DG 

805.27 

159.35 

697.17 

1129.62 

153.09 

178.15 

-DC 
“min  f2 

867.08 

155.61 

729.81 

1278.94 

147.66 

179.45 

-DC 

“min 

0.50 

868.61 

171.54 

641.68 

1095.52 

159.43 

183.64 

-DG 

936.58 

166.67 

701.72 

1171.47 

154.67 

178.53 

-DC 

“min/2 

1131.64 

162.99 

843.53 

1419.79 

150.45 

175.58 

-DC 
“min  A 

0.25 

1033.65 

172.95 

723.19 

1137.18 

156.55 

178.42 

-DG 

1076.53 

171.25 

743.61 

1187.43 

156.32 

176.24 

-DG 
“min  f2 

1207.33 

169.07 

835.23 

1331.34 

158.64 

173.47 

-DC 
“min  /] 

0.00 

1144.36 

179.03 

744.71 

1144.31 

163.82 

179.03 

-DC 

1197.79 

176.62 

749.21 

1197.74 

160.93 

176.62 

-DC 
“min  f2 

1307.33 

172.87 

828.55 

1307.35 

159.78 

172.87 

A 
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Fig.  12.  Bubble  plots  EENS  v/s  ECg.  Diameter  of  bubbles  proportional  to  CVaR(ENS)  (A)  and  CVaR(Cg)  (B). 


lower  for  increasing  /3.  This  to  be  expected:  lowering  the  values  of  / 3 , 
the  MOO  tends  to  search  for  optimal  allocations  and  sizing  ,EDC  that 
sacrifice  expected  performance  at  the  benefit  of  decreasing  the  level  of 
risk  (CVaR).  These  insights  can  serve  the  decision  making  process  on 
the  integration  of  renewable  DC  into  the  network,  looking  not  only  at 
the  give-and-take  between  the  values  of  EENS  and,  but  also  at  the  level 
of  risk  of  not  achieving  such  expected  performances  due  to  the  high 
variability. 

Fig.  13  shows  the  average  total  DG  power  allocated  in  the 
distribution  network  and  its  breakdown  by  type  of  DG  technology 
for  the  optimal  £°G  as  a  function  of  / 3 .  It  can  be  pointed  out  that  the 
contribution  of  EV  is  practically  negligible  if  compared  with  the  other 
technologies.  This  is  due  to  the  fact  that  the  probability  that  the  EV  is 
in  a  discharging  state  is  much  lower  than  that  of  being  in  the  other 
two  possible  operating  states,  charging  and  disconnected  (see  Fig.  10), 
combined  with  the  fact  that  when  EV  is  charging  the  effects  are 
opposite  to  those  desired. 

The  analysis  of  the  results  for  different  /3  values  also  allows 
highlighting  the  impact  that  each  type  of  renewable  DG  technology 
has  on  the  network  performance.  As  can  be  noticed  in  Fig.  13(A), 
the  average  total  renewable  DG  power  optimally  allocated,  increases 
progressively  for  increasing  values  of  /3:  this  could  mean  that  to 


obtain  less  'risky'  expected  performances  less  renewable  DG  power 
needs  to  be  installed.  However,  focusing  on  the  individual  fractions 
of  average  power  allocated  by  PV,  W  and  ST  (Fig.  13(B),  (C)  and  (E), 
respectively),  show  that  a  reduction  of  the  risk  in  the  EENS  and  ECg 
is  achieved  specifically  diminishing  the  proportion  of  PV  power 
(from  0.29^1  to  0.11^=0)  while  increasing  the  W  and  ST  (from 
0.38/J=1  to  0.48^=o  and  from  0.31/S=1  to  0.39^=o,  respectively),  but 
this  increment  of  W  and  ST  power  is  not  enough  to  balance  the  loss 
of  PV  power  due  to  the  limits  imposed  by  the  constraints  in  the 
number  of  each  DG  technology  to  be  installed  given  by  Tj.  Thus,  PV 
power  supply  is  shown  to  most  contribute  to  the  achievement  of 
optimal  expected  performances,  but  with  higher  levels  of  risk.  On 
the  other  hand,  privileging  the  integration  of  W  and  ST  power 
supply  provides  more  balanced  optimal  solutions  in  terms  of 
expectations  and  of  achieving  these  expectations. 

Table  7  summarizes  the  minimum,  average  and  maximum 
total  renewable  DG  power  allocated  per  node.  The  tendency  is  to 
install  more  localized  sources  (mainly  nodes  4  and  8)  of  renew¬ 
able  DG  power  when  the  MOO  searches  only  for  the  optimal 
expected  performances  (/?=1)  and  to  have  a  more  uniformly 
allocation  of  the  power  when  searches  for  minimizing  merely  the 
CVaR  (P=0). 
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Fig.  13.  Average  total  DG  power  allocated  (A)  and  its  breakdown  by  type  of  DG:  PV  (B),  W  (C),  EV  (D)  and  ST  (E). 


Table  7 

Average,  minimum  and  maximum  total  DG  power  allocated  per  node. 


PT  [kW]  P 


Node 

1.00 

0.75 

0.50 

0.25 

0.00 

Min 

Mean 

Max 

Min 

Mean 

Max 

Min 

Mean 

Max 

Min 

Mean 

Max 

Min 

Mean 

Max 

1 

12.08 

34.44 

54.77 

1.15 

22.40 

38.56 

0.00 

19.23 

40.98 

0.00 

39.03 

121.00 

3.00 

17.33 

34.71 

2 

2.30 

40.72 

69.73 

0.00 

49.95 

77.70 

36.50 

58.40 

123.36 

3.00 

63.61 

132.93 

0.00 

42.54 

84.09 

3 

0.00 

24.83 

46.45 

14.80 

41.79 

85.03 

0.00 

37.94 

105.11 

4.00 

36.87 

98.53 

1.00 

32.84 

77.78 

4 

76.00 

110.00 

133.41 

1.15 

67.40 

133.63 

0.58 

38.04 

80.13 

6.15 

20.73 

61.85 

0.00 

39.85 

85.86 

5 

22.60 

52.39 

77.08 

28.90 

60.66 

98.59 

12.63 

89.39 

143.50 

3.30 

23.49 

54.25 

1.00 

24.97 

79.64 

6 

12.33 

55.56 

85.46 

10.45 

21.22 

38.95 

2.00 

27.68 

106.26 

12.15 

53.78 

84.43 

0.00 

50.64 

116.85 

7 

8.00 

16.52 

35.38 

39.38 

64.07 

104.05 

0.00 

52.03 

159.73 

0.00 

34.09 

92.81 

5.00 

18.51 

39.23 

8 

79.03 

111.20 

146.63 

30.00 

74.57 

114.41 

0.00 

40.60 

146.06 

4.00 

37.94 

102.60 

1.00 

39.49 

119.38 

9 

0.00 

20.03 

68.73 

4.00 

74.07 

107.88 

0.00 

46.72 

85.61 

0.00 

44.06 

94.08 

0.00 

32.86 

74.53 

10 

0.00 

9.07 

25.35 

0.00 

1.58 

7.88 

0.00 

11.88 

58.69 

0.00 

8.58 

43.40 

0.00 

30.12 

83.45 

11 

0.00 

9.98 

17.68 

0.00 

3.04 

13.20 

0.00 

4.74 

23.45 

0.00 

8.99 

45.95 

0.00 

7.31 

51.17 

5.  Conclusions 

We  have  presented  a  risk-based  simulation  and  multi-objective 
optimization  framework  for  the  integration  of  renewable  generation 
into  a  distribution  network.  The  inherent  uncertain  behavior  of 


renewable  energy  sources  and  variability  in  the  loads  are  taken  into 
account,  as  well  as  the  possibility  of  failures  of  network  components. 
For  managing  the  risk  of  not  achieving  expected  performances  due  to 
the  multiple  sources  of  uncertainty,  the  conditional  value-at-risk  is 
introduced  in  the  objective  functions,  weighed  by  a  /?  parameter 
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which  allows  trading  off  the  level  of  risk.  The  proposed  framework 
integrates  the  Non-dominated  Sorting  Genetic  Algorithm  II  as  a 
search  engine,  Monte  Carlo  simulation  to  randomly  generate  realiza¬ 
tions  of  the  uncertain  operational  scenarios  and  Optimal  Power  Flow 
to  model  the  electrical  distribution  network  flows.  The  optimization 
is  done  to  simultaneously  minimize  the  energy  not  supplied  and 
global  cost,  combined  with  their  respective  conditional  value-at-risk 
values  in  an  amount  controlled  by  (3. 

To  exemplify  the  proposed  framework,  a  case  study  has  been 
analyzed  derived  from  the  IEEE  13  nodes  test  feeder.  The  results 
obtained  show  the  capability  of  the  framework  to  identify  Pareto 
optimal  sets  of  renewable  DC  units  allocations.  Integrating  the 
conditional  value-at-risk  into  the  framework  and  performing 
optimizations  for  different  values  of  j3  has  shown  the  possibility 
of  optimizing  expected  performances  while  controlling  the  uncer¬ 
tainty  in  its  achievement.  The  contribution  of  each  type  of  renew¬ 
able  DC  technology  can  also  be  analyzed,  indicating  which  is  more 
suitable  for  specific  preferences  of  the  decision  makers. 
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